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The zero temperature phase diagram of a one-dimensional S = 2 Heisenberg ferromagnet with 
single-ion cubic anisotropy is studied numerically using the density-matrix renormalization group 
method. Evidence is found that although the model does not involve quadrupolar couplings, there 
is a purely quadrupolar phase for large values of the anisotropy. The phase transition between the 
magnetic and quadrupolar phases is continuous and it seems to be characterized by Ising critical 
exponents. 
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I. INTRODUCTION 

It is well-known that crystal fields present in real mate-, 
rials can influence their magnetic behavior considerably.!!! 
In this paper we study a ferromagnetic insulator, with 
localized interacting magnetic moments, which are also 
subject to a crystal field of cubic symmetry. The 
finite-temperature critical phenomena in these systems 
have been investigated jpreviously within the Landau- 
Ginzburg- Wilson theoryjj and the cubic universality class 
and its critical exponents have been analyzed. However, 
in such phenomenological theories, the form of the order 
parameter is simply postulated, using mainly symmetry 
arguments. To see whether such an order really develops 
in the low-temperature regime one must start with the 
microscopic quantum Hamiltonian. 

In terms of a quantum spin Hamiltonian, the magnetic 
moments are described by spin operators, and the crystal 
field takes the form of a single-site term. The lowest value 
of S for which a cubic anisotropy term has any effect is 
5 = 2. Using this value our model Hamiltonian takes the 
form 



H 



(1) 



with Sf (a = x,y, z) denoting spin operators for S — 2 
at lattice site i. The first term is the isotropic, ferro- 
magnetic exchange interaction between nearest neigh- 
bors, and the second term represents the crystal field 
anisotropy; for 5 = 2 this is the most general single-site 
operator with cubic symmetry. The exchange term alone 
produces ferromagnetic order in the low-temperature 
phase; the question we want to address is how this order 
is affected by the cubic crystal field. 

The above Hamiltonian has been treated within j-the 
mean-field approximation (MFA) by many authors.!3 It 
has been found that the sign of the crystal-field con- 
stant D defines the easy axes of the magnetic ordering at 
low temperatures. For D = the system is isotropic and 



the spontaneous magnetization can lie along any direc- 
tion; for D < the cube diagonals [111] are preferred; for 
< D < (4z)/3, where z is the number of nearest neigh- 
bors, the directions parallel with the cube edges [100] are 
chosen. Moreover, the MFA predicts that for D > {4:z)/3 
the system is disordered, irrespective of the temperature. 

Recently, two oLus (M.D. and J.S.) have investigated 
the above modeljj applying the Raleyigh-Schrodinger 
perturbation theory in the limit D —^ 00. This calcula- 
tion is valid for large D, and the following picture emerges 
from it. In agreement with the MFA, for large D, there is 
no magnetic order at any temperature, i.e., the averages 
of the spin operators vanish: 



(5f)=0, for 



x,y,z. 



However, in contrast to the MFA prediction, the symme- 
try may still be spontaneously broken by the appearance 
of a purely quadrupolar order, described by the averages 

asm ^ asm = hsd')- 

(One may interchange x, y, and z in the above rela- 
tion.) This purely quadrupolar phase was predicted to 
exist in the model on a two-dimensionaUattice at least 
at zero temperature. In three dimcnsiongj and in higher 
dimensions it should be present both in the ground state 
and at finite (but low) temperatures. We suppose, that 
the quadrupolar phase also exists in one dimension at 
zero temperature (note that the broken symmetry is 
not continuous but discrete), although in this case the 
perturbation-theory arguments are weaker. 

The result that there may appear a purely quadrupo- 
lar order in the model in Eq. (|l|), which contains only the 
bilinear, "magnetic" exchange and the single-site term, is 
rather surprising. Usually, higher or|der exchanges are re- 
sponsible for quadrupolar ordering.!3 The appearance of 
quadrupolar order in the present case is entirely a quan- 
tum effect; one would not obtain purely quadrupolar or- 
der by simply replacing the spin operators in Eq. (|l|) 
with classical vectors, and then seeking the configuration 
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with minimum energy. The result is also of methodolog- 
ical interest. It is easy to understand that one cannot 
get purely quadrupolar order in this model by the MFA: 
only magnetization can act as a mean field when the ex- 
change term contains only first powers of the spin oper- 
ators. Therefore, if the purely quadrupolar phase exists, 
we have the example that, in contrast with the general 
belief, the MFA may give an answer which is qualitatively 
wrong in high dimensions. 

Unfortunately, there is no fully reliable method which 
would enable us to investigate the model in two or three 
dimensions, for general D. Therefore, the aim of the 
present paper is to check and confirm the existence of 
the purely quadrupolar phase numerically at least in the 
one-dimensional situation at zero temperature. We use 
White's xlensity-matrix renormalization group (DMRG) 
method,Ll which has a proven record of extreme efficiency 
in dealing with such problems. The method allows us to 
study the model for general D >Q, not only in the limit 
of large D as in Ref. ^, and we obtain the zero tempera- 
ture phase diagram, together with some properties of the 
(quantum) phase transition separating the magnetic and 
quadrupolar phases. 

The paper is organized as follows: we begin with sum- 
marizing the information obtained from the MFA and the 
perturbation theory in Sec. ^ Then we discuss how long- 
range order can be detected using the DMRG in Sec. Ell. 
Section IV is devoted to the concrete numerical results, 
while Sec. V is a brief discussion of our conclusions. We 
also attach an Appendix in which numerical predictions, 
obtained using identical techniques, are confronted with 
exact results in the test case of the ID Ising model in a 
transverse field. 



II. PREDICTIONS FROM MFA 
AND PERTURBATION THEORY 

Let us first list the eigenvalues and eigenstates jV'/c) 
of the single-site term in the Hamiltonian in Eq. ffl) 



ei = 


-241? 




62 = 


-24D 


1^2 


63 = 


-ISL* 


1^3 


64 = 




1^4 


65 = 




1^5 



(2) 



where by |2), . . . , | — 2) we denote eigenstates of . 
The five states are split into a doublet and a triplet, and 
the doublet energy is lower for Z? > 0. The doublet is 
nonmagnetic, i.e.. 



=0, 



(3) 



the nonmagnetic states are preferred in the total wave 
function of the system. 

For D — Q the model is isotropic and possesses the 
classical, ferromagnetic ground state with saturated mag- 
netization. As mentioned in the Introduction, the MFA 
predictaJ that the magnetic phase exists for < Z? < 8/3 
(z = 2 for the chain) , and that the directions of the cube 
edges [100] are the easy axes in this phase. According to 
the MFA, the ground-state magnetization decreases with 
D and disappears continuously at D = 8/3. The MFA 
result that the magnetic order vanishes at some finite D 
is a strong argument that this is indeed the case, since 
the MFA usually overestimates the tendency towards or- 
dering. Nevertheless, the actual value of D at which the 
magnetic phase ends may be much smaller than 8/3 (for 
example, in the Ising chain in a transverse field the MFA 
critical field is two times larger than the exact critical 
field) . Whether the actual phase transition is continuous 
or discontinuous is also a matter of question. We be- 
lieve that the easy axes in the magnetic phase are given 
correctly by the MFA. This approximation neglects the 
quantum fluctuations in the ground state, but it seems 
very unlikely that these could change the easy-axis di- 
rections, especially in the vicinity of D = 0, where the 
magnetization is almost saturated and the quantum fluc- 
tuations are small. 

Considering the large- 13 limit now, we make use of the 
perturbation theory developed in Ref. IJ. Since the two 
nonmagnetic states \'ipi) and |?/'2) dominate for large 
the model reduces to an effective two-state Hamiltonian. 
Restricting our attention to the one-dimensional situa- 
tion from now on, the effective Hamiltonian HcS can be 
written in terms of Pauli matrices af acting on the states 
\'ipi) and IV'2), and has the forma 

= --^ E (^^'^^i + -J^H ('^^^i) 

i i 



+(^(;^), (4) 



with 



~ 1 1 1 /I 

^ ^ 2D ' ~my^ ^ T2D^ ^ [Ij^ ) ' 



1 1 



o 



8D^ 321)3 XD-^ J ' 



for any i,j = 1,2 and a = x,y,z. Clearly, the crystal 
fleld with D > Q opposes the magnetic ordering, since 



(The exchange coupling constant J from Ref. ^ is set 
to 1 in the present work.) The quadrupole operators 
[{SfY — 2] of the original model in Eq. (^ turn out to 
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be-,represented by linear combinations of Pauli matrices 
aM 

[{5f)2-2] ^ ^/3af-< + 0(l/^?2), 

[{S:f-2] K-> 2a^ + 0{l/D'). (5) 

The analogous representations of 5"^ and {S" S'^ + 3'^' Sf), 
where a ^ /3, vanish up to 0{1/D^). 

In the lowest order in 1/-D, only Ji is present in 
the effective Hamiltonian, and Hcs describes the pla- 
nar model with continuous SO (2) symmetry in the XZ 
plane. While this model, conventionally written as the 
XY model, is ordered at zero temperature in two di- 
mensions and above, in one dimension it is known to be 
critical due to enhanced quantum fluctuations u'lI The in- 
plane correlation function 

decays asymptotically as a power law 1/Vl, where I is 
the distance between the two sites. Hence, on the basis 
of Eq. , we argue that for 13 cxd a critical state, in 
which the quadrupolar correlation function 

decays asymptotically as l/Vl, is approached. 

Taking into account the higher-order terms in HcS, 
we obtain the planar, XY-type model with small per- 
turbations. The three-spin term with the coupling con- 
stant J4 is crucial, because it reduces the symmetry of 
the effective Hamiltonian from the full rotational sym- 
metry in the XZ plane to the symmetry of discrete ro- 
tations through angles ±(27r)/3. In the one-dimensional 
case the Je-,and J3 terms are expected to be marginally 
irrelevant ,E3 but in the lack of a rigorous study of the off- 
diagonal three-point symmetry-breaking term with J4, 
we can only speculate (and then eventually check numer- 
ically) whether the 0{1/D^) contributions are relevant 
or irrelevant. If they are relevant, they can stabilize the 
"magnetic order" (i.e., the nonvanishing averages of af) 
in the effective model Hcs, since the Mermin- Wagner and 
Coleman's theorems do not hold when the symmetry is 
not continuous but discrete. If this is indeed the case, 
then Eq. ^ implies that for large D a quadrupolar or- 
der, associated with the operators JYS'f )^— 2] , exists in the 
original model of Eq. (|^) (see Ref. |4|for details). Alterna- 
tively, if all the perturbations turn out to be irrelevant, 
there should be a critical phase for large D. 

To summarize the above, we present the relevant (or- 
der) parameters to be analyzed as suggested by the MFA 
and the perturbation theory: 

m" = (Sr), = ((5^)2-2). (6) 

These ground-state averages are bulk values in the ther- 
modynamic limit. Note, that in a disordered state m" = 
= for a = x,y, z. Based on the above analysis we 



expect that for < D < Dm (where is probably 
smaller than the MFA value 8/3) a magnetic phase with 

771^ ^ 0, my = m' = 0, > = q\ (7) 

exists, while for D > _D,„ a purely quadrupolar phase 
with 

= my = m"- =0, q"" > qy ^ (8) 

emerges. (The indices x, y, z in the above descriptions 
of the phases may be interchanged.) In both phases the 
original cubic symmetry of the Hamiltonian is sponta- 
neously broken. Note however that the subgroup re- 
maining invariant is different in the two cases. As dis- 
cussed above, for large D a critical (disordered) phase 
with q" — cannot be excluded either. In any case, the 
values of must go to zero for D ^ 00. 

III. THE METHOD 

The density-matrix renormalization group (DMRG) 
method is one of the most reliable numerical techniques 
to study one-dimensional quantum lattice problems. For 
a detailed description of the algorithm see Ref.0. There 
are different implementations of the technique; neverthe- 
less it seems that the most accurate results can be ob- 
tained by investigating finite systems with open bound- 
ary conditions. In this setup the DMRG provides good 
approximations for the ground state and low-lying ex- 
cited states of long, but finite quantum chains. This can 
then be supplemented by a detailed finite size scaling 
analysis using appropriate scaling assumptions. 

Usually for the case of spontaneous symmetry break- 
ing the existence of long-range order is checked by study- 
ing appropriate two-point correlation functions {An An'). 
The direct observation of one-point functions {An) is 
apparently hindered by the fact that in finite systems, 
where tunneling is always finite between different vacua, 
the ground state is a symmetric combination of all pos- 
sible ordering directions, and thus naively {An) = 0. A 
possible way out, however, is to break the symmetry ar- 
tificially by a (small) auxiliary field which then forces the 
system to make a choice between the a priori undeter- 
mined directions, without noticably changing the value 
of the order parameter. This method has considerable 
advantages in the DMRG numerics where the measure- 
ment of two-point functions is rather awkward due to the 
open boundary condition. 

In our implementation the symmetry-breaking auxil- 
iary fields are only applied to the first and last spins of 
the open chain. These two boundary fields are chosen to 
be identical, leaving the system symmetric for i L — i 
reflections (L is the chain length). Then the Hamiltonian 
which is actually simulated in the numerical calculations 
is 

H ^ H ~h- (Bi+Bl) (9) 
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where and h are vectors of appropriate single-site op- 
erators and boundary fields. Since the system is not 
translation invariant, the one-point function (An) de- 
pends on the position n near the chain ends, but ap- 
proaches its bulk value in the middle of the system. If the 
system size is much larger than the correlation length the 
one-point function rapidly saturates and forms a plateau. 
Alternatively, if the chain length is too small comparing 
to the correlation length (e.g., the system is critical or it 
is close to criticality) no plateau appears to develop. 

There are two ways to estimate the order parame- 
ter and the associated correlation length from one-point 
functions. One can define the value 0^/2 = rnea- 
sured in the middle of the chain and study its depen- 
dence on the total chain length L. It is expected that 
o-L /2 (L) a (with a the bulk value of the order param- 
eter) exponentially fast as L — > 00, whenever the system 
is away from criticality, while the convergence is only al- 
gebraic at criticality. A disadvantage of this method is 
that many independent DMRG runs must be done with 
different L values, or one is forced to use the "infinite 
lattice" algorithm which is less precise. 

Alternatively, one can analyze the profiles an = (^n) 
as a function of n at a fixed length L, providing that L 
exceeds the correlation length ^ reasonably (which can 
be checked a posteriori) . In this latter case the profile is 
expected to have the following asymptotic form far from 
the chain ends 



g-(L+l-n)/^ 

'(L + 1 - n)x' 



(10) 



where a is the bulk value of the order parameter, ^ is the 
associated correlation length, x is a suitable exponent, 
and c is a constant. In general one can assume that the 
order parameter and the correlation length measured this 
way are identical to their bulk values derived from two- 
point correlation functions. However, the exponent of 
the polynomial prefactor x is normally different from its 
bulk value and may depend on the boundary conditions 
chosen. 

In practice we analyze our numerical data by dividing 
the chain into smaller segments with a length of the order 
of ^, and perform local least-squares fits using the formula 
in Eq. (10). We have found it appropriate to fix x first 



and then look for its value which makes the other fitting 
parameters the most stable in different segments as rt — > 
00. Our procedure of determining order parameters and 
correlation lengths was extensively tested on a simpler 
problem, the ID Ising model in a transverse field (ITF), 
where exact answers to many questions of interest are 
available. We summarize our experience with this model 
in the Appendix. 

In most cases the DMRG is very efficient to calculate 
spectral gaps. The method is less useful, however, if the 
ground state is (asymptotically) highly degenerate, since 
then several states have to be targeted together, resulting 
in a rapid drop in accuracy. In the cubic model under 
investigation we expect a 6-fold ground-state degeneracy 



in the magnetic phase, and a 3-fold degeneracy in the 
quadrupolar phase, in accordance with the number of 
easy-axis directions for the two order parameters. The 6- 
fold degenerate case already exceeded our computational 
limitations, thus some of our conclusions with respect to 
that case are based on short chain exact diagonalization 
results only. 

Errors in the quantitative predictions we report in the 
next sections stem from two sources, either from the lim- 
ited accuracy of the numerical technique we have been 
using, or from the approximate (asymptotic) nature of 
the formulas we applied for fitting and extrapolating the 
data. As for the DMRG errors, our results are some- 
what more precise when the "finite lattice" algorithm 
was used. In this case several iterations were done until 
convergence was reached. In principle, results obtained 
through the "infinite lattice" method inevitably contain 
a systematic "environment" error, and thus should be 
treated with less confidence. However, in the present 
problem improvements due to the "finite lattice" itera- 
tions turned out to be rather small, providing enough 
ground to assume that even our "infinite lattice" results 
have sufficient accuracy. 

In general we made various runs keeping different 
number of states AI, and extrapolated our results for 
M 00, or at least checked that the results had been 
converged in this parameter. Typically, we found good 
convergence in M despite the fact that the maximum 
number of states kept did not exceed M = 85. Note 
that the number of degrees of freedom at a single site is 
five in the S = 2 case, and the model in question has no 
continuous axial symmetry (total is not conserved) 
which could have been utilized to facilitate the DMRG 
calculation in the usual way. 



IV. NUMERICAL RESULTS 

In order to obtain the zero temperature phase diagram 
and critical properties of the model we carried out exten- 
sive numerical calculations using exact diagonalization 
techniques on short chains and the DMRG method on 
longer systems. In the case of open boundary conditions 
with boundary fields the actual form of the Hamiltonian 
used in the simulations was 

L-l L 



-h (Bi + Bl) 



(11) 



with a reduced cubic crystal field at the chain ends (to 
counterbalance the missing bonds for the first and the 
last spins) 



D, = 



D Hi 
D/2 ifi 



2,... 



1 



(12) 



and a general symmetry-breaking boundary field contain- 
ing magnetic and quadrupolar terms 
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/ h \ 
f u \ 


















hyy 







(13) 



In order to avoid handling complex numbers, we did not 
include a boundary term proportional to S"^. 

After the ground state had been found, we measured 
the order parameters. To facilitate further discussion we 
introduce the notations 



As was discussed in the former section, we could only 
apply the DMRG for the 3-fold degenerate case. We 
calculated the low-lying four states in an open chain at 
D = 1.3, which is above the critical point Dm (see later 
sections) , by targeting four states without any boundary 
fields. The results, presented in Fig. |l|(b), strongly sup- 
port the 3-fold asymptotic degeneracy with a finite gap 
A sa 0.08±0.02 above it in the thermodynamic limit. Un- 
fortunately, this calculation was rather time consuming, 
which impeded us to repeat it for other D values. 



< = (^">, q:^{{S:r-2) (14) 
with n = 1, . . . , L; a = x, y, z. 

A. Ground state degeneracies and gaps 

Figure |^(a) shows the lowest energy levels in a peri- 
odic chain with L — 8 sites obtained using exact diago- 
nalization techniques. Although the chain is short, the 
6-fold degeneracy, consisted of a singlet, a doublet, and a 
triplet is clearly discernible for small values of D. In the 
large D (quadrupolar) region the degeneracy seems to 
be 3-fold, consisted of the singlet and the doublet only. 
In both regimes we observe signs of a gap above these 
states, although the number of data points (L = 4,6,8) 
available for a given D did not allow us to perform a de- 
tailed finite-size scaling study. Note that since the bro- 
ken symmetry is discrete we do not expect any massless 
Goldstone modes in the ordered regimes. Excited states 
are rather massive domain walls, which separate domains 
with different ordering directions, and propagate in the 
system. 




FIG. 1. (a) Low-energy spectrum vs cubic anisotropy D 
in a short chain with L = 8 using periodic boundary con- 
dition. Shading schematically indicates higher-energy states, 
(b) Ground state-first excited state (singlet-doublet) gap (1st 
XS), and ground state-second excited state gap (2nd XS) 
as a function of 1/L (solid lines) and (dashed lines) 

at D = 1.3. The ticks in the horizontal axes correspond to 
L = 10, 20,..., 50. 



B. Magnetic order 

Let us now investigate the spontaneous magnetization 
and the corresponding correlation length. In these calcu- 
lations, purely magnetic boundary fields were used, i.e., 
hxx = hyy ^ hzz = in Eq. (|l|). 




L (b) 



FIG. 2. (a) The parameter m2/2 "^s the chain length L, 
and (b) the profile mj; vs the position n at L fixed, for some 
representative values of D. Purely magnetic boundary fields 
applied along the x axis as given by h^. 

Expecting the magnetic easy axes along the cube edges 
we started with applying a boundary field along the x 
axis setting hx — 10 or 0.1, ~ 0. Due to the symme- 
try, with such a boundary field one must obtain ^ 
and = — 0. Figure |^(a) shows how the param- 
eter to2/2i measured in the center of the chain, changes 
with the total chain length L for representative values 
of the crystal- field strength D. One can see that for 
D — 0.5,1,1.2 the parameter m^^^ converges to finite 
values independent of hx- This clearly indicates a finite 
magnetic order for these D values. On the other hand, 
for D = 1.3,1.5,2 the parameter m^/a S^es to zero as 
L —^ oo, meaning no magnetic order. As characteris- 
tic examples the magnetization profiles mjj with L fixed 
are shown in Fig. ^(b) for D = 1 and 2. For D = 1 a 
plateau corresponding to a finite spontaneous magneti- 
zation w 1.76, while for D = 2 a plateau with = 
can be observed. 

So far we have assumed that the cube edges give the 
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easy axes in the magnetic phase. To verify this as- 
sumption, we also apphed boundary fields which do not 
coincide with the expected easy directions. Choosing 
hx = 0.04 and hz — 0.03 we found that for L — > oo 
the parameter m^^^ always tends to the same value as 
found above, while to|^/2 go^^ to zero, in agreement with 



Dn 



1.2374 ±0.0004, 



/? = 0.127 ±0.004, 



Eq. (0). 




FIG. 3. The spontaneous magnetization m (left axis, dots) 
and the corresponding correlation length ^rn (right axis, tri- 
angles) as functions of D. Notice the large error bars for 
in the magnetic phase. 

The above calculations allowed us to determine how 
the bulk spontaneous magnetization m changes with the 
crystal- field strength D. We assume that monotonic de- 
pendencies of m.2/2 ^^^'^ those in Fig. ^(a), provide 
correct bounds for m, as was found for the ITF chain in 
the Appendix. The magnetization m is plotted in Fig. ^ 
as a function of D, together with the corresponding cor- 
relation length discussed in details below. Since £,m 
increases around the phase transition, very long chains 
have to be studied in order to find the large-L limit of 
m|^^2 i'^ this region. The largest D for which we can 
see the magnetic order is D = 1.23. For this D, tak- 
ing L = 400, we obtain a finite m with an uncertainty 
of 0.003 (for smaller D, shorter chains are sufficient to 
calculate m much more precisely). Due to the increasing 
correlation length, and thus computational limitations, 
we were not able to trace the dependence of m on Z? 
further. Starting from D = 1.26, the phase with no mag- 
netic order is observed. 

Although it cannot be shown directly that the spon- 
taneous magnetization disappears continuously at the 
phase-transition point, the sharp increase of the corre- 
lation length suggests that the transition is continuous. 
We therefore attempt to fit the dependence of m on D, 
shown in Fig. ^, assuming a power law singularity and a 
linear regular term 

m{D) = ci{Drn - Df + C2{D^ - D). (15) 

The fit is very good and we obtain 



ci = 2.15 ±0.05, C2 = -0.17 ±0.05. Interestingly, the 
exponent /3 is very close to that of the two-dimensional 
Ising model. The obtained parameter £).,„ is our best 
estimate of the phase-transition point. 

To calculate the magnetic correlation length 5,„, we 
analyze the profiles of at fixed L, analogous to those 
presented in Fig. §(b). As was described earlier, we di- 
vide the chain into short sections and perform local fits 
using the formula in Eq. (10). 



For each D > Dm, i.e., in the phase with no magnetic 
order, the local values of the fitting parameters obtained 
with X = are very stable, and the convergence with 
increasing n (towards the center of the chain) is convinc- 
ing. The fitted order parameter m is always very close 
to zero. It seems that the decay of is asymptoti- 
cally purely exponential, like the decay of rn}^^ in the 
disordered phase of the ITF chain, discussed in the Ap- 
pendix. Due to the stability of local ^'s, the correlation 
length 5m in the phase with no magnetic order can be 
calculated very precisely (at D — 1.27 the error in 5m is 
around 0.1% and it is smaller for larger D). 

In the magnetically ordered phase, for D < D„i, the 
local fitting parameters obtained with x = are not sta- 
ble. The local 5's increase towards the center of the chain, 
like in the ordered phase of the ITF chain, suggesting the 
existence of an algebraic prefactor. To calculate 5m, we 
vary x and look for its value where the local 5's are most 
stable. We do not find a universal prefactor, and the ob- 
tained X depends on D and h^- However, the estimates 
of 5m are independent of the boundary field, which jus- 
tifies the procedure; the uncertainties in 5m are around 
10%. (The same procedure was tested for the ordered 
phase of the ITF chain and yielded acceptable estimates 
of the correlation length — see Appendix.) 

The dependence of the correlation length 5m on the 
crystal-field strength D is plotted in Fig. ^. Due to the 
large error bars for D < Dm, the functional dependence 
can only be analyzed for D > Dm- Here, assuming a 
second order transition again, we fit the dependence of 
5m on D with a power law 

^m{D)^C3{D~D'mr''. (16) 

The fit is again rather good, yielding the values 

Dm = 1.236 ±0.004, = 1.02 ±0.06. 

Thus, close to the phase transition both the spontaneous 
magnetization and the corresponding correlation length 
seem to be well described by power laws. The two inde- 
pendent estimates of the phase-transition point, Dm and 
D'm from Eqs. (|l^) and (|l^), nicely coincide with each 
other. (In what follows, we will use Dm, which is more 
precise.) Moreover, the exponents /3 and i' are close to 
those of the two-dimensional Ising model. Altogether, we 
have found strong indications that, as far as the magnetic 
quantities m and 5m are concerned, the phase transition 
is continuous and is described by the critical exponents 
of the two-dimensional Ising model. 
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C. Quadrupolar order 

We now investigate the quadrupolar order, associated 
with the parameter defined in Eq. (^). While the 
quadrupolar order must be present in the magnetic phase 
(D < Dm), the question whether it exists in the phase 
with no magnetic order {D > Dm) is of central interest. 

Expecting the quadrupolar easy axes along the cube 
edges, as given in Eqs. (0) and (||), we began with apply- 
ing boundary fields along the x axis. These were purely 
quadrupolar: ha = 0; h^^ 7^ 0, hyy = hzz = in Eq. (13); 
or purely magnetic: hx ^ 0, hz = 0; haa = 0. With such 
boundary fields, due to the symmetry, one must obtain 

<fu^ll = Qn- 
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FIG. 4. (a) The parameter 52/2 the chain length L, and 
(b) the profile qf^ vs the position n at L fixed, for two values 
of D. Purely magnetic (hx) and purely quadrupolar {hxx) 
boundary fields applied along the x axis. 

We found that in both phases the order parameter ^2/2 
measured in the chain center converges with increasing 
L to finite values which do not depend on the actual 
boundary fields. Figure ^(a) shows these dependencies 
for D = 1.2 < Dm and D = 1.3 > Dm- The appropriate 
profiles of with L fixed are plotted in Fig. ^(b) . In both 
phases the picture is typical for an ordered state: has 
plateaus in the middle of the chain, which correspond 
to finite bulk q". Thus, we find quadrupolar order on 
both sides of Dm, confirming the existence of a purely 
quadrupolar phase for D > Dm- 

To verify that the cube edges give the easy axes of the 
quadrupolar ordering, we applied a quadrupolar bound- 
ary field having all three components different from zero: 
ha = 0, = 0.03, hyy = 0.02, and h^z = 0.01. Then, 
in both phases, it was observed that g£^2 > 1^/2 ~ 1l/2 
in the large-L limit, and the limiting values of 92/2 '^^^^ 
the same as found above, in accordance with Eqs. (0) 
and (|). 

Assuming that monotonic dependencies of g£^2 011 ^ 
provide correct bounds for the bulk value q, we deter- 
mined q for various values of D. (q is defined as the 
largest of q" in a broken-symmetry state.) The results 



are depicted as bold points in Fig. || the precision in q is 
10"'^ or better. We observed that, starting from D « 1.3, 
the associated correlation length increases rapidly with 
D (see the inset in Fig. |^). At the same time, the DMRG 
precision for q" fell dramatically. For these reasons, al- 
ready at L* = 1.45 (and for larger D too) we were unable 
to obtain precise results for a chain long enough to show 
us the large-L limit of q£^2 (for D = 1.45 and 1.5 we 
obtained upper bounds on q, which are shown as open 
circles in Fig. ||). We had similar problems close to the 
transition point Dm- 




FIG. 5. The quadrupolar order parameter g as a function 
of D. For D = 1.45 and 1.5 upper bounds on q are depicted 
as the open circles. The inset shows very rough estimates of 
the corresponding correlation length ^q. 

The effective Hamiltonian Hos in Eq. (^ may help to 
understand the difficulties for Z) > 1.4. Notice, that the 
symmetry-breaking term J4 in Hcs is very small, J4/ Ji = 
0{1/D^). Thus even if J4 is relevant, the correlation 
length is expected to diverge rather rapidly as I? — > 00. 
On the other hand, the observed abrupt increase in 
and the emerging numerical difficulties beyond Z? ^ 1.4 
may indicate an additional phase transition in which the 
quadrupolar long-range order vanishes, and gives rise to 
an extended critical phase with algebraically decaying 
quadrupolar correlations beyond a critical value Dq. A 
naive linear fit to the data points for q in Fig. |^ would 
yield a value Dq « 2.1, but we are not in a position 
to give any firm claims in this question. Note however 
that if Dq = 00, i.e., if the quadrupolar phase extends 
to Z? ^ 00, the curve q{D) should finally bend upwards, 
and in the data points in Fig. ||, together with the upper 
bounds at D — 1.45 and 1.5, the (small) curvature is 
consistently downwards. 

Now we analyze the dependence of g on Z) in Fig. ^ 
close to the phase-transition point Dm- For D < Dm we 
assume a power-law singularity and a linear regular term 



q{D) = qicft + CiiDn 
Fitting yields 



D) 



D) 



(17) 



gicft = 1.838 ± 0.01, 



= 0.5 ±0.1, 
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C4 = 0.44±0.1, C5 = -0.36±0.1. Thus, on this side of the 
phase-transition point, the quadrupolar order parameter 
appears to have singular behavior, but the exponent 6 
differs from [3 describing the spontaneous magnetization. 

For D > D„i, we suppose the behavior is regular, which 
is suggested by Fig. |5|. Our fitting formula is now 

q{D) = gright + c&{D - D^) + c^{D - D^f, (18) 

giving 

g.ight = 1.840 ± 0.002, 

C6 = -1.83 ± 0.03, C7 = -0.6 ± 0.1. The left and right 
estimates of the order parameter q at Dm are in agree- 
ment, which supports our scaling assumptions. However, 
we cannot exclude the possibility that a singular deriva- 
tive in q exists also for D > D,n, so close to that we 
do not observe it in our numerics. 

To calculate the correlation length ^g, we analyze the 
profiles of g", like those in Fig. ^(b). In our opinion, only 
a purely quadrupolar boundary field is suitable for this 
purpose. With a magnetic boundary field, the decay of 
could be driven by the decay of m", and so governed 
by and not ^q. 

Using the formula in Eq. again, we find that both 
for D < D„i and D > D,n the prefactor exponent x 
is nonzero and it depends on D and the boundary field. 
The fitting yields rough estimates to £,q (Fig. || inset), the 
errors are believed to be around 20%. The results suggest 
that — for D < which is the expected answer, 
since the magnetic and quadrupolar order is intimately 
connected in this phase. For D > Dm, the correlation 
length increases with D starting from D « 1.3, as 
mentioned above, but, due to the poor precision, we can- 
not resolve whether diverges (numerically it increases 
sharply) for D — s- Dm or it remains finite. 

V. CONCLUSIONS 

In summary, we have studied the ground-state prop- 
erties of a one-dimensional 5* = 2 ferromagnetic spin 
chain with single-site cubic crystal-field anisotropy D. 
We argued that in contrast with the mean-field predic- 
tion, perturbation theory suggests the possibility that a 
purely quadrupolar phase exists for large values of D. 
This quadrupolar phase was expected to be separated 
from the magnetic phase by a continuous quantum phase 
transition at a critical point Dm- We verified this con- 
jecture by investigating the model numerically using the 
density-matrix renormalization group method on open 
chains with special, symmetry-breaking boundary condi- 
tions. 

In most cases, the method provided precise estimates 
of the magnetic and quadrupolar order parameters. Very 
close to the phase-transition point and in the large-Z? 
limit, however, our results were less accurate due to the 
rapidly increasing correlation lengths. For the correlation 



lengths themselves we could only obtain rather rough es- 
timates. 

Evidence has been obtained that, in qualitative agree- 
ment with the mean-field prediction, the spontaneous 
magnetization diminishes continuously at the phase- 
transition point. Regarding the magnetic properties, the 
transition seems to be characterized by the critical expo- 
nents of the two-dimensional Ising model. This could be 
plausible, since the extra broken symmetry of the mag- 
netic ground state with respect to that of the quadrupolar 
ground state is just a Z2 subgroup — the same group as 
the one breaking down spontaneously in the case of the 
Ising model. 

Both in the magnetic and nonmagnetic phases we 
demonstrated the presence of a quadrupolar order. In 
the former case the quadrupolar order is just a "sec- 
ondary" effect, inevitably present in any S > 1 model 
with magnetic order. In the latter case, however, it is 
the "primary" order, constituting a qualitatively different 
broken-symmetry phase. At the phase-transition point 
Dm, the quadrupolar order parameter is continuous. For 
D ^ Dm, we could clearly discern a singularity in the 
derivative of the order parameter, and found that the 
quadrupolar correlation length diverges together with the 
magnetic correlation length. For D > Dm, our observa- 
tions are much less concrete: the order parameter can be 
fitted reasonably well by a low-order (regular) polynomial 
without any singular terms, and, although the quadrupo- 
lar correlation length increases as Dm is approached from 
above, our results do not suggest unambiguously a di- 
vergence on this side. It is not clear for us whether the 
phase transition involving the magnetization should have 
any precursor in the quadrupolar fiuctuations above Dm- 

Mainly due to computational limitations we were un- 
able to resolve convincingly the question whether the 
quadrupolar phase extends to Z) ^ cx3 or there is a finite 
value Dq where the quadrupolar long-range order disap- 
pears. This question is unique to the one-dimensional sit- 
uation — in two dimensions and higher we expect a finite 
quadrupolar order at zero temperature for any D > Dm- 

Finally, the scenario for the magnetic-to-quadrupolar 
phase transition was supported by the analysis of the de- 
generacy of the ground state as a function of D. We ob- 
served a 6-fold and a 3-fold asymptotic degeneracy below 
and above Dm respectively, in accordance with the the 
number of possible ordering directions in the two phases. 

The confirmation of the purely quadrupolar phase in 
the one-dimensional case gives rise to the belief that such 
a phase should also exist in higher dimensions at suf- 
ficiently low temperatures when the cubic crystal field 
is strong. Increasing temperature necessarily destroys 
the quadrupolar order and leads to a finite-temperature 
phase transition between the quadrupolar and a com- 
pletely disordered phase. Note that the correct critical 
theory of this transition cannot be obtained by simply 
substituting spins with classical vectors in our model, 
since this approach is unable to account for the purely 
quadrupolar order. The usual practice of neglecting 
quantum fluctuations by treating spins classically around 
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a finite-temperature phase transition would confront fun- 
damental difficulties in this case. 
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APPENDIX: A TEST CALCULATION 

Our numerical procedure and extrapolation methods 
were tested on the zero-temperature ID Ising model in 
a transverse field (ITF) where, an exact solution for the 
bulk quantities is available .EiliJ Our experience with this 
model provided useful guidelines to refine our procedure 
in the study of the cubic model of Eq. (^ . 

The Hamiltonian of the ITF model on an open chain 
with length L is defined as 



ITF 



L-1 



(Al) 



i=l 1=1 

F if i = 2,...,L-l 

F/2 if i = 1,L 



where af denotes Pauli matrices at site i, and F is 
the transverse magnetic field. Note that in the above 
definition the spins at the chain ends, being coupled 
to one neighbor only, are subject to the reduced field 
F/2. We concentrate on the spontaneous magnetization 
j^iTF _ ^^2^ ^-^^ -j-j^g corresponding correlation length 
^ITF^ For F = the Hamiltonian H^^^ has the classi- 
cal, ferromagnetic ground state with mF^^ = il, and the 
single-site term in H^"^^ tends to destroy the spontaneous 
magnetization; in this respect the ITF model is similar to 
the cubic model in Eq. (|l|). In the thermodynamic limit 
there are two phases in the ITF chain, distinguished by 
the order parameter mF^^ . For < F < 1 the ground 
state is ordered due to spontaneous symmetry breaking 



± (l - F 



r-ITF 



1 



21ogF' 

while for F > 1 there is a disordered phase with 

1 



^ITF ^ 0^ ^ITF 



logF 



(A2) 



(A3) 



ITF 



At F = 1 the ground state is critical and ^ 

In our study the following values of the transverse field 
were considered: F = 0.9 (ordered phase) where S}^^ ~ 
4.75; F = 1.1 (disordered phase) where S}^'^ w 10.49; and 
the critical point F = 1 with ^^"^^ = oo. We computed 



the position dependent order parameter mj^ 
using the DMRG method with boundary fields up to 
L = 200. As a symmetry-breaking boundary field we 
considered 



B = a^ 



(A4) 



m Eq. (|), and we used two different values of ft-: a strong 
/i = 10 and a weak ft = 0.1 boundary field. In the case 
of the ITF chain the DMRG is extremely precise even 
close to the critical point .E3 Keeping a relatively small 
number of states M — the truncation errors are al- 
ready negligible and the fitting procedure can be tested 
on numerically exact data. We emphasize that although 
rigorous results are available for quantities in the ther- 
modynamic limit, the behavior of relatively short chains 
with complicated boundary conditions is not feasible for 
a study with purely analytical tools even for the ITF 
model. 
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FIG. 6. ITF chain results, (a) Order parameter 
sured in the center of the chain vs the chain length L. (b) 
Order parameter profiles function of position in the 

chain n at fixed L = 200. Curves are labeled by the values of 
the transverse field F and the boundary-field strength h. 



The value of the order parameter at one of the central 
spins ^ ^ function of the chain length L is shown 

in Fig. ^(a). For these data the "infinite-lattice" DMRG 
was used from L = 4 to L = 200, thus a systematic 
"environment error" cannot be excluded.E3 In the ordered 
phase, converges to the same finite value for both 

values of ft., and the two dependencies are monotonic . A t 
each L, correct bounds for the bulk mP-^ in Eq. ( [A^ ) 
are obtained, and at L = 200, ten digits of m}^^^ are 
recovered. In the disordered phase, ti^j converges to 
zero for both boundary fields ft. At the critical point, 
the limit of does not show up in Fig. ||(a), we can 

only say that the correlation length must be of the order 
of L = 200 or larger. 

In Fig. ^(b) we show the order parameter profiles m}-^^ 
for a chain length fixed at L = 200. This length exceeds 
the exact correlation lengths for F = 0.9 and F = 
1.1 many times. There are plateaus for F = 0.9 and 
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r = 1.1, which correspond to the bulk values of the order 
parameter. There is no plateau at the critical point where 
the correlation length is infinite and the dependence of 
the profile is algebraic. 

In order to calculate the correlation length S}"^^ , we 
analyze the profiles from Fig. ^b) using the ansatz in 
Eq. (10). As is described in Sec. Ill we first fix a value 
of the exponent x ^^id then fit for the other parameters. 
Repeating the fits in different sections of the chain we 
seek the value of x which makes the other parameters 
the most stable as a function of n. 

For r = 1.1 the value of the exponent which yields 
the most stable fitting parameters is x « 0. From po- 
sitions n = 50 to n = 100, and for both h, the fitted 
correlation length agrees with the exact £}'^^ with a 
precision of 7 digits! It seems that, asymptotically, the 
decay of m^n^ in the disordered phase is purely expo- 
nential. This should be compared with the exponent 
Xbuik ~ 1/2 characterizing the decay Cj£ the bulk two- 
point correlation function in this phase.tl (For a similar 
difference in algebraic prefactors of one and two-point 
correlation functions in an S' = 1 chain see Ref. O.) 



Anticipating computational limitations in the investi- 
gation of the 5 = 2 cubic model in Eq. (m), we also 
analized shorter chains, such that L « 10^^""" during the 
test calculations. For F = 1.1 and L = 100 four digits 
of the exact ^^"""^ are recovered; for F = 0.9 and L = 50 
only rough estimates for ^i^^^ g^j.g obtained, the error is 
around 10 %. 

Our observations can be summarized as follows: 

• The order parameter m}"^^ can be calculated accu- 
rately in both phases; it can be decided whether 
the system is ordered or not. 

• In the ordered phase, upper and lower bounds for 
m}^^ are obtained by considering strong and weak 
boundary fields, respectively. 

• The correlation length ^^'^^ can be calculated very 
precisely in the disordered phase. 

• In the ordered phase, rather rough estimates of ^^'^^ 
can be found. 




n 



FIG. 7. ITF chain results. Ordered phase F = 0.9. The 
local values of the correlation length ^ from Eq. (p^), ob- 
tained with different boundary fields h and exponents as a 
function of the position n at L = 200. The exact correlation 
length ^^^^ is shown as a horizontal line. 

In the ordered phase for F = 0.9 we can also start with 
X = 0. In this case, however, the fitting parameters are 
not stable, and ^ increases with n for both values of h, 
as is illustrated by Fig. ^ This suggests that a nonzero 
prefactor is present. The value of x which makes the 
correlation length ^ the most stable seems to depend on 
the boundary field applied. We obtain x ~ 1-25 ± 0.05 
for h — 10 and x ~ 0.9 ± 0.05 for h = 0.1, we cannot see 
a universal prefactor. (Note that for the bulk two-point 
functions the exponent of the prefactor is Xbuik = 2 in 
this phase.) The values of f calculated with the above 
prefactors are shown in Fig. 0. It is seen that using these 
prefactors we find better estimates for ^^"""^ than with 

x = o. 
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